Self-organized criticality in the hysteresis of the Sherrington - Kirkpatrick model 
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We study hysteretic phenomena in random ferromagnets. We argue that the angle dependent 
magnetostatic (dipolar) terms introduce frustration and long range interactions in these systems. 
This makes it plausible that the Sherrington - Kirkpatrick model may be able to capture some 
of the relevant physics of these systems. We use scaling arguments, replica calculations and large 
scale numerical simulations to characterize the hysteresis of the zero temperature SK model. By 
constructing the distribution functions of the avalanche sizes, magnetization jumps and local fields, 
we conclude that the system exhibits self- organized criticality everywhere on the hysteresis loop. 



PACS numbers; 64.60.Lx, 75.60.Ej, 75.10.Nr 

Hysteresis in ferromagnetic systems is a century old 
physical problem. Efficient phenomcnologies have al- 
ready been developed |l|], but an accepted microscopic 
theory is yet to be constructed. In soft magnets, where 
domain wall motion dominates the physics, considerable 
progress has been achieved recently [^-Q. In hard mag- 
nets domain nucleation, domain wall motion and their 
interaction are all important. Hence they are better de- 
scribed on a more microsopic level as an assembly of 
strongly interacting spins or hysterons |q|. Quantitative 
insight to such systems has been gained recently through 
studying the random field Ising model (RFIM) H . 

However a key aspect of the physics of real systems is 
missing from the RFIM: it does not include the long range 
dipolar (or magnetostatic) interactions. While these are 
negligible on atomic scales relative to the exchange term, 
they can dominate the collective behaviour of granular 
systems. This is so because the dipolar interaction is 
long ranged, so it involves every spin in the volume of the 
grains, whereas the exchange coupling scales only with 
the number of spins on the surface of the grain. These 
dipolar forces are important: they prevent the rough- 
ening of the domain walls B and determine the size of 
the domains Q| . Crucially, the sign of these interactions 
changes with the angle. This introduces frustration into 
the system, which is not represented in the RFIM. 

To capture the influence of frustration on hysteretic 
phenomena, we study the simplest system, contain- 
ing long-range frustrated interactions, the Sherrington- 
Kirkpatrick (SK) model. Early numerical work demon- 
strated that this model exhibits hysteresis [0j|]. How- 
ever, in spite of its obvious importance, we could not find 
analytic studies of the hysteresis loop of the SK model. 
In this Letter we use scaling arguments, replica calcula- 
tions and large scale numerical simulations to character- 
ize the hysteresis of the zero temperature SK model. By 
constructing the distribution functions of the avalanche 
sizes, magnetization jumps and local fields, we conclude 
that the system exhibits self-organized criticality every- 
where on the hysteresis loop. 



The SK model consists of N Ising spins (a = ±1) on a 
fully connected lattice, described by the Hamiltonian 
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where Jij is a random Gaussian number of zero mean and 
variance 1/iV. Throughout the paper we work at T = 0. 

First we summarize our numerical results. We start 
from a fully polarized state and change the external mag- 
netic field h adiabatically: for a given field we let all spins 
align according to their local field before varying h again. 
During the avalanches we use sequential single spin fiip 
updating to ensure the decrease of the total energy. The 
resulting hysteresis loop for the SK model is presented 
in Fig. [I]. Finite size scaling analysis shows that the hys- 
teretic trajectories are well-defined in the A^ — > cx) limit, 
and the coercive field converges to a finite value. 

We also analyzed the minor hysteresis loops of the SK 
model (inset Fig. |lj). Within numerical accuracy they re- 
turn to the major loop at the point of departure, exhibit- 
ing return point memory. This feature is present in many 
experimental systems, and it is also one of the criteria for 
the applicability of the Preisach phenomenology [|l| . 

Next we establish some of the basic energy scales from 
elementary considerations. When spin aj is flipped, the 
local field hi at another site changes by an amount pro- 
portional to 2Jij ^ 2/N^/^. Thus the external field h 
has to be changed by an amount dh ex 1/N^'^ to start 
a new avalanche. Now let S be the change in the total 
magnetization during an avalanche, and dm — S/N the 
jump of the magnetization rn during the avalanche. The 
average m{h) curve is continuous and thus its deriva- 
tive (dm/dh) ex {S)/N^/'^ is finite (Fig. |^), requiring: 
(5) ex N^/^. This is possible only if the scale of the 
distribution of avalanches is set by iV^'^. This is char- 
acteristic of systems at criticality, whereas off-criticality 
the scale is set by some control parameter of the Hamil- 
tonian. This leads to the central result of the paper: the 
SK model exhibits critical behaviour everywhere along 
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FIG. 1. The hysteresis loop of the SK model, averaged 
over 100 disorder configurations (A'^ = 1600). Inset: multiple 
minor loops, exhibiting return point memory. 



its hysteresis loop. As this phenomenon is independent 
of the parameters of the Hamiltonian, it is a manifesta- 
tion of self-organized criticality. 

To elucidate this point, in Fig. we show the distri- 
bution functions of S, and the number of spin flips in an 
avalanche (its "size"), n; V{S) and ^{n), respectively, 
measured in the interval m g [—0.3,0.3] for various sys- 
tem sizes. Both distributions exhibit power law behavior 
and can be well described by the finite size scaling forms: 



V{n) = (B/lnN) n'^ din/N") 
V{S) = {A/\nN) S-^ p{S/nP) 



(2) 
(3) 



withr, e= 1±0.1, cr = 0.9±0.1,and/3 = 0.6±0.1. The 
logarithmic prefactors were necessary to achieve satisfac- 
tory scaling collapse. Since such terms are needed only 
to keep distributions with an exponent 1 normalized, this 
strongly suggests that t = g — \ exactly. Unfortunately, 
because the cutoffs of the distributions 7^(5') and I?(n) 
scale with different powers of N the attractive picture of 
a diffusive motion of the local fields due to the random- 



ness in Ji. 



would lead to an infinite diffusion constant 



D ex (n)/N^/^ and is thus inapplicable. 

Adopting the r = 1 equality and combining it with 
(S) ~ N^'^ immediately yields the relation (3 — 1/2, 
with logarithmic corrections, in good agreement with the 
above measured value. Also, because the Jij's take neg- 
ative values as well, spins of both signs are destabilized 
in an avalanche. Therefore the number of participating 
spins is only bounded from below by S'/2, yielding the ex- 
ponent - bound a > P = 1/2. An upper bound for a can 
be obtained from estimating the dissipated energy, Ed, 
during a finite but small sweep of the external-field hi -^ 
h2 = hi + Ah: Ed = NmAh ~ N. Also, since the aver- 
age energy dissipation per spin is at least 2dh ^ 1/iV^/^, 
Ed can be estimated as Ed > 2ntotai/N^/'^, where ritotai 
is the number of flips during all avalanches from hi to 
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FIG. 2. Avalanche size and magnetization jump distribu- 
tions V{n) and V{S) for system sizes N = 400, 800, 1600, 
and 3200. The inserts show the collapsing scaling curves 
corresponding to Eq. (H) with d{n/N'^) — nhi{N)'D{n) and 
p{S/N'^) ^Sln{N)P{S). 



/i2. But the number of avalanches during this sweep is 
proportional to Ah/dh ^ N'^/'^, i.e. ntotai ^ N^/^(n). 
Combining all the above gives N ^ Ed > (n) ^^ N"' im- 
plying the upper bound cr < 1 , which is nearly saturated 
according to our numerics. 

The above distributions imply that the average value of 
X = dm/ dh is dominated by a few very large avalanches, 
whereas its typical value scales to zero as ~ 1/iV^/^, which 
we confirmed independently numerically. Therefore the 
hysteresis loop for a specific disorder realization has a 
slope zero with unit probability, interrupted by a few 
macroscopically large avalanches. This feature is char- 
acteristic of the Barkhausen noise and establishes the 
frustrated spin glasses as possible candidates to describe 
certain classes of hysteretic magnets. 

We also studied the correlations of consecutive 
avalanches. We measured the Hausdorf dimensions of 
the numerically determined hysteresis loop and that of a 
sequence of independent avalanches, generated with the 
above distributions. Having found the two Hausdorf di- 
mensions equal suggests that avalanches are uncorrelated. 

These results, in particular the size A'' as the sole cutoff 
of the different distribution functions, which all exhibit 
power law behaviour, confirm the above - stated self- 
organized criticality of the entire hysteresis loop of the 
SK model. To shed more light on the underlying physics 
we explore the local-fields, hi — ^ Jij<^j + ^ by study- 




FIG. 3. The distribution of the local stabilities, P(A), for 
iV = 3200. Inset: The finite size scahng of P(0). 



ing the local stabilities, Ai = aihi, which are all positive 
for stable spin configurations. Their distribution, P{X) is 
shown in Fig. |^. Remarkably - unlike the local field dis- 
tribution Q ~ P{X) is essentially the same at any point 
of the hysteresis loop. This suggests that the avalanche 
dynamics of the SK model organizes the system into spe- 
cial states with similar properties everywhere along the 
hystersis loop. A careful finite size analysis shows that 
P{X = 0) - 1/VN and P(A) « CX" with C « a « 1 for 
small A's. As we now show, this latter result establishes 
once again that these special states are critical. To prove 
this let us flip rifup arbitrary spins starting from a given 
stable spin configuration {fXi} with Ai > and calcu- 
late the average number of new unstable spins, {riunst): 
distinguished by negative stabilities A^ = A^ + AA^ < 0: 



A,- — A,; 



^i*Ji' 



j flipped 



(4) 



The system is critical if (n„„st) = n/Hp, as for (uunst) < 
rifiip the avalanches die out exponentially fast while in 
the opposite case they explode Q. Assuming that the 
n/iip random terms at the rhs. of Eq. (j^ are indepen- 
dent, the probability Pd of destabilizing a given spin is: 



Pd = 



dX P{X) / d{AX) g(AA) 



(5) 



where Q{AX) = cxp{—NAX'^/8nfUp}y/N/87rnfii„ is the 
probability distribution of the A A term in Eq. (Q), and 
P(A) is approximated by its asymptotic form, P(A) = 
CX". The average number of destabilized spins is then 
(riunst) = NPd = C(a)iV(n/,,p/iV)("+i)/2, with Cia) an 
a-dependent constant, C(l) = C. For a > 1 (or a = 1 
and C < 1) (riunst) < nfup and the system cannot give 
rise to large avalanches. On the other hand, for a < 1 
(or a — 1 and C > 1) (riunst) > nfUp, and the state is 
unstable. Thus the criticality condition is characterized 
by a = 1 and C = \. These are exactly the values found 
in our numerical simulations, once again underlining the 
criticality of the system. 



The physical mechanism of self-organized criticality 
can be qualitatively understood as follows. As the 
avalanche rolls, at any given time step t the stabilities of 
the spins are shifted only by those spins, which changed 
sign at step t — 1. These spins have flipped because the 
second term of Eq. (^) for their stabilities was negative, 
pulling their A^'s downward. However once Xi changed 
sign, the very same term now enhances this stability. 
More importantly, this term being symmetric, it also 
pushes upward the stabilities of the other spins of the 
avalanche, which pulled spin i down and flipped it in the 
first place. This effect is suppressing the density of states 
with low local fields, reminescent of the formation of the 
Coulomb gap in the disordered electron problem |9[] . The 
stabilities of the spins which did not participate in the 
avalanche will be shifted by a random amount by the 
just-flipped spins. However in the presence of a slope in 
their distribution P{X), this will have a net effect, mov- 
ing the stabilities of more spins downward than upward. 
In short, correlations between the spins of an avalanche 
move the stabilities of the already flipped spins upward; 
at the same time the random couplings between all spins 
drive a net downward drift. The competition of these two 
forces keeps the system critical. 

To understand the shape of the measured major hys- 
teresis loop more in detail we first observe that the states 
where an avalanche stops must always be single spin-flip 
stable (SSS). Let us therefore define the average number 
of SSS states. 
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(V{m,h)) = lT^{\{Q{Xi) 5{mN -Y^a,)] 



(6) 



where the angular bracket indicates annealed disorder 
averaging and the trace stands for the summation over all 
spin configurations. The product of the theta functions 
in Eq. (0) projects out those states where all the spins 
have positive stabilities, while the delta function selects 
states with a given magnetization. 

Using the integral representations 5{y) — /^ ^^-ixy 
and Q{Xi) = JZc^^e-^'-^' the function V{m,h) 
can be rewritten in an exponential form, V{m, h) = 
I ^UtU ^ exp{-i(A,Zi - xa.^ + mx)}] , and the dis- 
order average and the spin summation can be easily car- 
ried out. After the disorder averaging the effective action 
contains a term proportional to ~ (Si-^*)^- Decoupling 
this term with a new Hubbard-Stratonovich field R, one 
finally arrives at the following expression: 

(^) = ^/^/ ^ / ^ exp{iV[lng - :^ - zmx]} , 

Q(/i, x,p)= / ^ 2 cos(a; - zh) e'^''-'^' . (7) 

J 27rz 

The above integral can readily be evaluated in the N -^ 
(X) limit by the saddle point method. The saddle point 
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FIG. 4. Outer bound of the region of the single spin-flip 
stable states. Dotted line: Jo = 0; solid line: Jo = 2.5. The 
impossiblity of a monotonic m{h) curve within these bounds 
forces the jump, as indicated by the arrows. 



equation, dQ/dx = imQ can be solved analytically, and 
the variable x can be completely eliminated resulting in 
the following expression for {V{m,h)): 
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{V{h,m)) r^exp{Nnsp{m,h)) , 
1 - m . l-d)- 1 
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where R is determined from d ftsp{ni, h, R)/d R — 0, and 
^± = 0(^) with <j>{x) = -^ /; e-^'dt. 

In Fig. 4. we plotted the contour of flgp = 0. Out- 
side this line the density of SSS states scales to ex- 
ponentially, thus they are definitely unable to arrest the 
avalanches. Inside this line the number of SSS states is 
exponentially large, and is thus comparable to the total 
number of states, themselves expoential in N. Therefore 
avalanches get trapped with a higher probability in one 



of the SSS states. Hence the fig 



contour constitutes 



a strict outer bound for the true hysteresis loop. 

Comparing Figs.l and 4 shows that the 51 ^p = con- 
tour considerably overestimates the size of the hysteresis 
loop. We pursued two refinements of this calculation. 
We developed a replica symmetric description as well as 
a de Almeida - Thouless type replicon instability analy- 
sis: these will be reported separately |10|] . 

Finally, we briefly discuss the effect of a finite ferro- 
magnetic coupling, Jo > 0. Jo simply shifts the value of 
the magnetic field h ^ h + Jom in Eq. (M), and results 



in a shear of the entire contour 17, 



0. For a branch 



of the hysteresis loop ?77, must be a monotonic function 
of h. Since the flgp — contour is an outer bound of the 
hysteresis loop, therefore, when this loop would force a 
non-monotonic m{h) relation (Fig. 4), the major hystere- 
sis loop must develop a finite jump. Since the slope of the 
major hysteresis loop for Jo = is finite, one expects this 
transition to occur at a finite critical coupling, Jq ~ Jc- 
Our numerical data agree with this picture pOl . 

We end with a comparison to the random field Ising 
model. In that model our initial simple scaling consid- 



erations yield (S) oc 0(1), i.e. a non-critical avalanche 
distribution. As shown in Ref. |^, there is only a single 
critical point at the coercive field at some specific value of 
the disorder. Therefore the distribution functions exhibit 
a scaling behaviour, with the cutoff set by the distance 
from this critical point, rather than by the system size, as 
happens for the SK model. In short, the RFIM exhibits 
"plain old criticality" M, while the SK model exhibits 
self-organized criticality. Also, the avalanche distribu- 
tion exponents r in the two models are different: on the 
mean field level r = 1.5 for the RFIM and 1.0 for the 
SK model. In finite dimensions the exponents typically 
increase: numerical simulations of the 3D RFIM found 
T = 1.6 [|6|. In contrast, numerical studies of realistic 
3D models |1^ , as well as experimental [Q works report 
T in the 1.1. ..1.4 regime. This raises the possibility that 
the finite dimensional extensions of the frustrated models 
might provide a r closer to the experimental values. 

In sum, we studied the hysteretic behaviour of the SK 
model. We determined numerically the distribution func- 
tions of the avalanches, the magnetization jumps, and the 
local fields. The model exhibits self-organized criticality 
everywhere along the hysteresis loop. We recalculated the 
loop with analytic methods as the location of one-spin- 
fiip stable states, and found satisfactory agreement with 
the numerical results. 
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